/*******************************************************************************

nbody-view.c copyright 2003 Joseph MacMillan

NBodyView reads in positions from an N-body data file and displays the
particles using OpenGL.

usage:  NBodyView [options] infile
        options:
        -W #        Initial window width
        -H #        Initial window height
        -i file     Use file for input
        -x #        Rotate by # degrees around x-axis
        -y #        Rotate by # degrees around y axis
        -n #        Only use first # of particles
        -m #        Mouse tolerance, in pixels.
        -o          Use ortho projection, rather
                    perspective

Controlling movement:
        <up>        Translate upward
        <down>      Translate down
        <left>      Translate to the left
        <right>     Translate to the right
        <pgup>      Zoom in (translate along z axis)
        <pgdown>    Zoom out (translate along z axis)
        <x>         Rotate by 1 degree around x axis
        <X>         Rotate by -1 degree around x axis
        <y>         Rotate by 1 degree around y axis
        <Y>         Rotate by -1 degree around y axis
        <z>         Rotate by 1 degree around z axis
        <Z>         Rotate by -1 degree around z axis
        <space>     Reset rotations and translations to default
        <q, Q, esc> Quit

        <left mouse button>     Rotate
        <Right mouse button>    Translate
        <Middle mouse button>   Zoom in/out

*******************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <getopt.h>
#include <GL/gl.h>
#include <GL/glext.h>
#include <GL/glu.h>
#include <GL/glut.h>

//
// Global variables
//

// Window information
int width = 600;            // width of window
int height = 600;           // height of window
int win;                    // GL window ID

// Projection to use
int ortho = 0;

// mouse info
int left_mouse_moving = 0;
int right_mouse_moving = 0;
int middle_mouse_moving = 0;
int mouse_start_x = 0;
int mouse_start_y = 0;
int mouse_tolerance = 1;

// GL tranformations
float rot_x = 0.0;         // rotation about x-axis
float rot_y = 0.0;         // rotation about y-axis
float rot_z = 0.0;         // rotation about z-axis

float trans_x = 0.0;       // translation about x-axis
float trans_y = 0.0;       // translation about y-axis
float trans_z = 0.0;       // translation about z-axis

// Use vertext extension only if direct rendering is available
int use_vertex = 0;

// N-Body system info
int num_parts;              // Total number of particles to plot
double snap_time;           // Time of snap shot
float *p;                   // pointer to array of positions

float r_max;                // maximum radius -- for drawing limits

// glut specific functions
// Plots the particle positions, using current rotations and tranlations.
void cbPlotPoints(void);
// Called when a window is dynamically resized.
void cbResizeWindow(int w, int h);
// Called when a ker is pressed.
void cbKeyPressed(unsigned char key, int x, int y);
// Called when a "special" key is pushed.
void cbSpecialKeyPressed(int key, int x, int y);
// Called when a mouse button is clicked
void cbMouseClicked(int button, int state, int x, int y);
// Called when the mouse moves with a button pressed down
void cbMouseMovement(int x, int y);

// Reads in data, checks command line args, etc
void Init(int argc, char *argv[]);
// Does a few gl initializations
void InitGL(void);
// Get command line arguments
void GetArgs(int argc, char *argv[], char *filename, int *test, int *n_ask);
// Prints command line help
void PrintOptions(void);
// Prints movement help
void PrintUsage(void);
// Specialized function to set up *p
void LoadParticles(FILE *stream, int n_ask);
// To test plotting
void LoadTestParticles(void);

int main(int argc, char *argv[])
{
    // I'm using GLUT here for ease of use.
    glutInit(&argc, argv);

    // Set up stuff
    Init(argc, argv);
    
    // Want rbga colours, float buffering, and depth
    glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);
   
    // Set initial window size
    glutInitWindowSize(width, height);

    // Create the window.  win is the window ID
    win = glutCreateWindow("NBodyView");

    // Call this is do the drawing
    glutDisplayFunc(&cbPlotPoints);

    // Call this when the window is resized
    glutReshapeFunc(&cbResizeWindow);

    // Handle key presses
    glutKeyboardFunc(&cbKeyPressed);
    glutSpecialFunc(&cbSpecialKeyPressed);
    
    // Handle mouse events
    glutMouseFunc(cbMouseClicked);
    glutMotionFunc(cbMouseMovement);
    
    // a few GL initializations
    InitGL();
    
    // Main loop -- the program never gets control back
    glutMainLoop();

   return 0;
}

void Init(int argc, char *argv[])
{
    int n_ask = -1;
    int test = 0;
    char filename[80] = "";
    FILE *fp;
    
    fprintf(stderr, "NBodyView: plot particle postions using OpenGL\n");      
    
    // Get command line arguments
    GetArgs(argc, argv, filename, &test, &n_ask);
   
    // Load in particles
    if (test)
        LoadTestParticles();
    else if (strlen(filename) == 0)
        LoadParticles(stdin, n_ask);
    else
    {
        if ((fp = fopen(filename, "r")) == NULL)
        {
            fprintf(stderr, "Can't open file %s for reading\n", filename);
            exit(EXIT_FAILURE);
        }
        LoadParticles(fp, n_ask);
        fclose(fp);
    }
    
    //PrintUsage();
    
}

void InitGL(void)
{   
    // set up initial placement of viewer
    trans_z = -3.0 * r_max;
        
    // The gl stuff
    glClearColor(0.0, 0.0, 0.0, 0.0);
        
#ifdef GL_EXT_vertex_array
    {
        const GLubyte *test_render = glGetString(GL_RENDERER);
            
        if (strstr(test_render, "Indirect") == NULL)
        {
            use_vertex = 1;
            glEnable(GL_VERTEX_ARRAY_EXT);
        }
    }
#endif
        
    // Set up some other gl stuff
    cbResizeWindow(width, height);
}

void GetArgs(int argc, char *argv[], char *filename, int *test, int *n_ask)
{
    int option;
    int NumArgs = 0;
    while ((option = getopt(argc, argv, "hi:W:H:n:x:y:m:ot")) != EOF)
    {
        switch (option)
        {
            case 'i':
                strncpy(filename, optarg, 80);
                NumArgs += 2;
                break;
            
            case 't':
                *test = 1;
                NumArgs += 1;
                break;

            case 'W':
                width = atoi(optarg);
                NumArgs += 2;
                break;
                
            case 'H':
                height = atoi(optarg);
                NumArgs += 2;
                break;
            
            case 'o':
                ortho = 1;
                NumArgs += 1;
                break;
            
            case 'n':
                *n_ask = atoi(optarg);
                NumArgs += 2;
                break;
            
            case 'x':
                rot_x = atof(optarg);
                NumArgs += 2;
                break;
                
            case 'y':
                rot_y = atof(optarg);
                NumArgs += 2;
                break;
            
            case 'm':
                mouse_tolerance = atoi(optarg);
                NumArgs += 2;
                break;
                
            case 'h':
            case '?':
                PrintOptions();
                PrintUsage();
                exit(EXIT_SUCCESS);
        }
    }
    
    if (argc > NumArgs + 1) strncpy(filename, argv[NumArgs + 1], 80);
}

void PrintOptions(void)
{
    fprintf(stderr, "NBodyView: use OpenGL to plot particles\n");
    fprintf(stderr, "Usage: NBodyView [options] infile\n");
    fprintf(stderr, "   [options]\n");
    fprintf(stderr, "   -W #        Initial window width\n");
    fprintf(stderr, "   -H #        Initial window height\n");
    fprintf(stderr, "   -i file     Use file for input\n");
    fprintf(stderr, "   -x #        Rotate by # degrees around x-axis\n");
    fprintf(stderr, "   -y #        Rotate by # degrees around y axis\n");
    fprintf(stderr, "   -n #        Only use first # of particles\n");
    fprintf(stderr, "   -m #        Mouse tolerance, in pixels\n");
    fprintf(stderr, "   -o          Use ortho projection, rather than");
    fprintf(stderr,               " perspective\n");
}

void PrintUsage(void)
{
    fprintf(stderr, "Controlling movement:\n");
    fprintf(stderr, "   <up>        Translate upward\n");
    fprintf(stderr, "   <down>      Translate down\n");
    fprintf(stderr, "   <left>      Translate to the left\n");
    fprintf(stderr, "   <right>     Translate to the right\n");
    fprintf(stderr, "   <pgup>      Zoom in (translate along z axis)\n");
    fprintf(stderr, "   <pgdown>    Zoom out (translate along z axis)\n");
    fprintf(stderr, "   <x>         Rotate by 1 degree around x axis\n");
    fprintf(stderr, "   <X>         Rotate by -1 degree around x axis\n");
    fprintf(stderr, "   <y>         Rotate by 1 degree around y axis\n");
    fprintf(stderr, "   <Y>         Rotate by -1 degree around y axis\n");
    fprintf(stderr, "   <z>         Rotate by 1 degree around z axis\n");
    fprintf(stderr, "   <Z>         Rotate by -1 degree around z axis\n");
    fprintf(stderr, "   <space>     Reset rotations and translations");
    fprintf(stderr,               " to default\n");
    fprintf(stderr, "   <q, Q, esc> Quit\n");
    fprintf(stderr,                "\n");
    fprintf(stderr, "   <left mouse button>     Rotate\n");
    fprintf(stderr, "   <Right mouse button>    Translate\n");
    fprintf(stderr, "   <Middle mouse button>   Zoom in/out\n");
}

void cbResizeWindow(int w, int h)
{
    // to avoid crash ?
    if (h == 0)
        h = 1;
    
    // set new window width and height
    width = w;
    height = h;
    
    // set viewing area (all of window)
    glViewport(0, 0, width,  height);
    
    // Set up perspective
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    
    if (ortho)
        glOrtho(-1.2 * r_max, 1.2 * r_max, -1.2 * r_max, 1.2 * r_max,
                 0.1 * r_max, 100.0 * r_max);
    else
        gluPerspective(45.0, (GLfloat)width / (GLfloat)height, 
                       0.1 * r_max, 100.0 * r_max);
    
    glMatrixMode(GL_MODELVIEW);
}

void cbPlotPoints(void)
{      
    // set colour to white -- can make this user selectable?
    glColor3f(1.0, 1.0, 1.0);
       
    // clear colour buffer
    glClear(GL_COLOR_BUFFER_BIT);
    
    // set matrix
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    // Transform to current
    glTranslatef(trans_x, trans_y, trans_z);
    glRotatef(rot_x, 1.0, 0.0, 0.0);
    glRotatef(rot_y, 0.0, 1.0, 0.0);
    glRotatef(rot_z, 0.0, 0.0, 1.0);
    
    // plot points
#ifdef GL_EXT_vertex_array

    if (use_vertex)
    {
        glVertexPointerEXT(3, GL_FLOAT, 0, num_parts, p);
        glDrawArraysEXT(GL_POINTS, 0, num_parts);
    }
    else 
#endif
    {
        int i;
        
        glBegin(GL_POINTS);
        for (i = 0; i < num_parts; i++)
        {
            glVertex3fv(&p[3 * i]);
        }
        glEnd();
    } 
    
    glFlush();
    
    glutSwapBuffers();
}
        
void cbKeyPressed(unsigned char key, int x, int y)
{
    switch (key)
    {
        case 113:  // q
        case 81:   // Q
        case 27:   // esc
            glutDestroyWindow(win);
            exit(EXIT_SUCCESS);
        break;
        
        case 'x':
            rot_x += 1.0;
            cbPlotPoints();
            break;
        case 'X':
            rot_x -= 1.0;
            cbPlotPoints();
            break;
            
        case 'y':
            rot_y += 1.0;
            cbPlotPoints();
            break;
        case 'Y':
            rot_y -= 1.0;
            cbPlotPoints();
            break;
               
        case 'z':
            rot_z += 1.0;
            cbPlotPoints();
            break;
       case 'Z':
            rot_z -= 1.0;
            cbPlotPoints();
            break;
            
        case 32:
            rot_x = rot_y = rot_z = 0.0;
            trans_x = trans_y = 0.0;
            trans_z = -3.0 * r_max;
            cbPlotPoints();
            break;
            
        default:
            fprintf(stderr, "No action for key %d.\n", key);
            break;
    }
}

void cbSpecialKeyPressed(int key, int x, int y)
{
    switch (key)
    {
        case GLUT_KEY_PAGE_UP:
            trans_z += 0.05 * r_max;
            cbPlotPoints();
            break;
            
        case GLUT_KEY_PAGE_DOWN:
            trans_z -= 0.05 * r_max;
            cbPlotPoints();
            break;
            
        case GLUT_KEY_UP:
            trans_y += 0.05 * r_max;
            cbPlotPoints();
            break;

        case GLUT_KEY_DOWN:
            trans_y -= 0.05 * r_max;
            cbPlotPoints();
            break;

        case GLUT_KEY_LEFT:
            trans_x -= 0.05 * r_max;
            cbPlotPoints();
            break;

        case GLUT_KEY_RIGHT:
            trans_x += 0.05 * r_max;
            cbPlotPoints();
            break;

        default:
            fprintf(stderr, "No action for %d.\n", key);
            break;
    }
}

void cbMouseClicked(int button, int state, int x, int y)
{
    if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN)
    {
        left_mouse_moving = 1;
        mouse_start_x = x;
        mouse_start_y = y;
    }
    if (button == GLUT_LEFT_BUTTON && state == GLUT_UP)
    {
        left_mouse_moving = 0;
    }
    
    if (button == GLUT_RIGHT_BUTTON && state == GLUT_DOWN)
    {
        right_mouse_moving = 1;
        mouse_start_x = x;
        mouse_start_y = y;
    }
    if (button == GLUT_RIGHT_BUTTON && state == GLUT_UP)
    {
        right_mouse_moving = 0;
    }
    
    if (button == GLUT_MIDDLE_BUTTON && state == GLUT_DOWN)
    {
        middle_mouse_moving = 1;
        mouse_start_x = x;
        mouse_start_y = y;
    }
    if (button == GLUT_MIDDLE_BUTTON && state == GLUT_UP)
    {
        middle_mouse_moving = 0;
    }
}

void cbMouseMovement(int x, int y)
{
    if (left_mouse_moving)
    {   
        rot_x += 180.0 * (y - mouse_start_y) / height;
        rot_y += 180.0 * (x - mouse_start_x) / width;
        
        if (abs(x - mouse_start_x) > mouse_tolerance || 
            abs(y - mouse_start_y) > mouse_tolerance)
        {
            cbPlotPoints();
            mouse_start_x = x;
            mouse_start_y = y;
        }
    }
    
    if (right_mouse_moving)
    {   
        trans_x += 2.0 * r_max * (x - mouse_start_x) / width;
        trans_y -= 2.0 * r_max * (y - mouse_start_y) / height;
        
        if (abs(x - mouse_start_x) > mouse_tolerance || 
            abs(y - mouse_start_y) > mouse_tolerance)
        {
            cbPlotPoints();
            mouse_start_x = x;
            mouse_start_y = y;
        }
    }
    
    if (middle_mouse_moving)
    {   
        trans_z -= 2.0 * r_max * (y - mouse_start_y) / height;
        
        if (abs(x - mouse_start_x) > mouse_tolerance || 
            abs(y - mouse_start_y) > mouse_tolerance)
        {
            cbPlotPoints();
            mouse_start_x = x;
            mouse_start_y = y;
        }
    }
}
        
void LoadTestParticles(void)
{
    int i, j;
    float radius;
    
    num_parts = 10000;
    
    // allocate memory
    p = malloc(3 * num_parts * sizeof(float));
    r_max = 0.0;
    for (i = 0; i < num_parts; i++)
    {
        for (j = 0; j < 3; j++)
        {
            p[3 * i + j] = -10.0 + 20.0 * rand() / (float)RAND_MAX;
        }
        radius = sqrt(p[3 * i + 0] * p[3 * i + 0] +
                      p[3 * i + 1] * p[3 * i + 1] +
                      p[3 * i + 2] * p[3 * i + 2]);
        r_max = (radius > r_max) ? radius : r_max;
    }
}

void LoadParticles(FILE *stream, int n_ask)
{
    int i, c, mod;
    float dump;
    float radius;
    float x, y, z;
    
    // first is number of particles -- integer
    if (fread(&num_parts, sizeof(int), 1, stream) != 1)
    {
        fprintf(stderr, "LoadParticles: can't get number of particles\n");
        exit(EXIT_FAILURE);
    }
    
    if (n_ask > 0 && n_ask < num_parts)
    {
        mod = (int)ceil((double)num_parts / (double)n_ask);
        n_ask = (int)ceil((double)num_parts / (double)mod);
    }
    else
        n_ask = num_parts;
    
    // second is time of snapshot -- may need this later
    if (fread(&snap_time, sizeof(double), 1, stream) != 1)
    {
        fprintf(stderr, "LoadParticles: can't get time\n");
        exit(EXIT_FAILURE);
    }
    
    // allocate enough memory to hold positions only
    if ((p = malloc(3 * n_ask * sizeof(float))) == NULL)
    {
        fprintf(stderr, "LoadParticles: can't get enough memory\n");
        exit(EXIT_FAILURE);
    }
    
    // read in just positions
    // warning! I'm doing away with error checking for this loop!
    mod = 1;
    if (n_ask < num_parts)
        mod = (int)ceil((double)num_parts / (double)n_ask);
    c = 0;
    for (i = 0; i < num_parts; i++)
    {        
        // first mass
        fread(&dump, sizeof(float), 1, stream);
        // next, three positions
        fread(&x, sizeof(float), 1, stream);
        fread(&y, sizeof(float), 1, stream);
        fread(&z, sizeof(float), 1, stream);
        // next, three velocities
        fread(&dump, sizeof(float), 1, stream);
        fread(&dump, sizeof(float), 1, stream);
        fread(&dump, sizeof(float), 1, stream);
        
        if ((i % mod) == 0)
        {
            p[3 * c + 0] = x;
            p[3 * c + 1] = y;
            p[3 * c + 2] = z;
            
            // do r_max calulation now
            radius = sqrt(p[3 * c + 0] * p[3 * c + 0] +
                          p[3 * c + 1] * p[3 * c + 1] +
                          p[3 * c + 2] * p[3 * c + 2]);
            r_max = (radius > r_max) ? radius : r_max;
            c++;
        }
        
        
   }
   num_parts = n_ask;
   fprintf(stderr, "Successfully read %d particles at time %f\n", 
           num_parts, snap_time);
}
